Phonon-affected steady-state transport through molecular quantum dots 
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We consider transport through a vibrating molecular quantum dot contacted to macroscopic leads 
acting as charge reservoirs. In the equilibrium and nonequilibrium regime, we study the formation of 
a polaron-like transient state at the quantum dot for all ratios of the dot-lead coupling to the energy 
of the local phonon mode. We show that the polaronic renormalization of the dot-lead coupling 
is a possible mechanism for negative differential conductance. Moreover, the effective dot level 
follows one of the lead chemical potentials to enhance resonant transport, causing novel features 
in the inelastic tunneling signal. In the linear response regime, we investigate the impact of the 
electron-phonon interaction on the thermoelectrical properties of the quantum dot device. 
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I. INTRODUCTION 

Electronic devices featuring a single organic molecule 
as the active element, so called molecular junctions, are 
promising candidates in the search for further miniatur- 
ization and novel functionality. Such systems can be de- 
scribed as quantum dots: mesoscopic systems coupled to 
macroscopic charge and heat reservoirs. 

Molecular junctions are susceptible to structural 
changes when being occupied by charge carriers. The 
local interaction with optical phonons becomes apparent 
as vibrational signatures in the current- voltage character- 
istics of the device [l|-[||, resulting from the interference 
of elastic and inelastic tunneling processes and the renor- 
malization of the effective dot level energy ■ When 
the vibrational energy and the electron-phonon (EP) in- 
teraction become sufficiently large, nonlinear phenomena 
emerge, such as hysteresis, switching and negative differ- 
ential conductance (NDC). As is well known from the 
Holstein molecular crystal model [H, , strong EP in- 
teraction may heavily reduce the "mobility" of electrons 
through the formation of small polarons himj. Thus, 
the formation of a local polaron is considered a possible 
mechanism for the observed nonlinear transport proper- 
ties of molecular junctions [lj| . 

Molecular junctions may also constitute efficient power 
generators or heat pumps. Their highly energy- 
dependent transmission together with the tunable level 
energy could be used to optimize the thermoelectrical 
figure of merit. In the weak dot-lead (DL) coupling 
limit, the theoretical efficiency approaches the Carnot 
value (l6| . However, long electron residence times in- 
crease the effective EP coupling. Moreover, some level 
broadening is needed to ensure useable power output. 



That is why, for practical applications, the regime of 
comparable electronic and phononic time scales becomes 
interesting. 

In our work, we calculate the steady-state charge and 
energy transport through the quantum dot for small- 
to-large DL coupling and weak-to-strong EP interac- 
tion. Based on a variational Lang-Firsov transforma- 
tion fl5l [T^-HH , we determine the nonequilibrium dot 
spectral function in the formalism of Kadanoff-Baym [22[ 
and calculate the dot self-energy in a self-consistent way 
up to second order in the renormalized interaction coeffi- 
cients. The variational parameter is determined numeri- 
cally by minimizing the thermodynamic potential. 



II. MODEL 

We consider the standard Hamiltonian of the single- 
site quantum dot. It is based on a modified Fano- 
Anderson model with the static impurity being replaced 
by a single site coupled to a local phonon mode (H = 1): 



H = (A - fi)dU - guorfdft + b)+ u tfb 



(1) 
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The quantum dot is represented by the energy level A, 
with the fermionic operators . It is coupled to a local 
phonon mode of energy u>o, with the dimensionless 
EP coupling strength g. The operators c£t? (for k = 
1, . . . , N; a = L,R) correspond to free electrons in the 
N states of the left and right lead, with the energies Ek a 
and the equilibrium chemical potential fj,. The last term 
in Eq. (JXJ) allows for dot-lead particle transfer. 

To account for the competition between polaron lo- 
calization and charge transport, we apply to the model 
an incomplete Lang-Firsov transformation [2l|, in- 
troducing the variational parameter 7 G [0,1]. Then 
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H = S\HS 7 , with 



5 7 = exp{< ? (6 t - b)(-y<fld + (l - 7 )n d )} 



(2) 



For 7 = 1, S-y coincides with the shift-transformation of 
the Lang-Firsov small-polaron theory fl7j . which elim- 
inates the EP coupling term in Eq. (JT]) and lowers the 
dot level by the polaron binding energy e p = g 2 coo- For 
7 < 1, S* 7 accounts for the quasistatic displacement of 
the equilibrium position of the oscillator, which is pro- 
portional to the dot mean occupation n d — (d'd). After 
the transformation the Hamiltonian reads 



H = fjtfd - gu (l-j)(b^ + b)(dU-n d ) 

+ LJ tfb + £p (l - jfnj + J2(tka - Ma C ka (3) 



In ([3]), the DL coupling is affected by the EP interaction. 
Furthermore the bare dot level is renormalized: 



Tf = A-fi- e P 7(2 - 7) - 2e p (l - 7) 2 n d 



(4) 



Note that now d and b are the operators of dressed elec- 
trons (in analogy to polarons) and the shifted oscillator. 
The original electron and oscillator operators now read 
d = exp{^g(tf — b)} d and b = b + jgd^d + (1 — j)gn d . 

The application of a potential difference between the 
leads is described by adding to (j3|) the interaction with 
the external fields U a = —Sfi a and defining the voltage 
bias with e being the negative elementary charge: 



int 



$ =(U L -U R )/e. (5) 



III. THEORETICAL APPROACH 

A. Polaronic spectral function in the 
Kadanoff-Baym formalism 

For finite voltage bias between the nonintcracting 
macroscopic leads, the response of the quantum dot is 
given by the polaronic noncquilibrium real-time Green 
functions 



9dd(tutr,U) = i(dl J (t 2 )d u (t 1 )) , 
g% d (h,h;U) = -iiduitjdlih)) 



where the time dependence of d^ 



(6) 
(7) 



l v is determined by 
H + Hi nt . According to Kadanoff-Baym [22j], the real- 
time response functions may be deduced using the equa- 
tions of motion for the nonequilibrium Green functions 
G dd (ti,t 2 ] U, t ) of the complex time variables t = t — iT, 
t € [0,/9]. We base our calculations on the Dyson equa- 
tion of the polaronic Green functions, which defines the 



polaronic self-energy Y, dd = G dd * — ^dd ■ F° r a gi yen 
ordering of t\, t 2l the equations of motion of the func- 
tions g dd follow through the limiting procedure to — > — oo. 
Limiting ourselves to the steady-state regime, we sup- 
pose that all functions depend only on t = t\ — t 2 . After 
a Fourier transformation by the method used in Ref. I22I 
the following exact equations for the steady-state are ob- 
tained HH: 

U)X> d (^ U) - <&(w; U)X< d (cj; U)=0, (8) 

[u}-rj-ReZ dd (u;U)]A{u;U) = 

r(w; U) Re g dd (uj; U) . (9) 



Here we defined, in analogy to the equilibrium case, 

A(u;U)=g> d (u;;U)+g< d (Lj;U), 
duj A{uj;U) 



gdd(z; U) 

T(oj;U) 
^dd(z; U) 



2ir z — ui 
r(u>-,U) = X> d (u,;U)+X< d (u;-,U) 

2tt z — uj 



(10) 

(11) 
(12) 

(13) 



where A(w; U) is the polaronic nonequilibrium spectral 
function. According to Eq. (|10[) . we can write 



9dd(^U) 
9dd(^ U) 



A(u; [/)/>; U) , 
A(u;U)(l-f(u;U)) 



(14) 
(15) 



introducing the noncquilibrium distribution /, which fol- 
lows from Eqs. ([8]) and (fl~2|) as 



/>;[/) 



(16) 



For the Green function g dd in Eq. (|lll) we use the ansatz 
gdd{z] U) = l/(z — rj — Yt dd (z; U)) and find the following 
formal solution of Eq. ©: 



A(lj;U) = 



T(u;U) 



-v-vf 



2k 



r(u>';U) 



(17) 



To deduce a functional differential equation for the self- 
energy T^dd we add to H mt in Eq.([5]) the interaction with 
fictitious external fields {V} (cf. Refs. fl~9l - [23h . The equa- 
tions of motion of the polaronic Green functions are then 
expressed by means of the functional derivatives of Tidd 
with respect to {V}. The resulting equations for EJ- rf are 
solved iteratively to the second order in the renormalized 
EP and DL interaction coefficients in ([3]), while the cor- 
relation functions of the interaction coefficients are eval- 
uated supposing independent Einstein oscillators. We 
then let {V} —> and perform the limit to —> — oo. A 
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subsequent Fourier transformation yields 
£<>;[/) = 4^ ( W ; (7) + [(1 - 7)5^ 



A(w - u ; U)f(u - uj q ; U)n B (u ) 
+A(lo + io a ; U)f(u + loo; U){n B {us Q ) + 1)1 , (18) 



,(i)<, 



(u;U) = J2 {-To(«)f (°)(w + ft)n F (u + U a ) 

a 

+ ^2Is(k)2 sinh(s0) IrjM (u - soj q + /1) 

S>1 

xn B (suj )n F (w - sw + f7 ) + r^ 0) (w + sw + fi) 

(19) 



x(ns(sw ) + l)n F (u + scjq + U a ) |, 
with = (e^ + l)" 1 , n B (w) = (e? u - 1 

2a M 



and 



e -vVcothe r (o) (a;)) 
2i:\t a {uj)\ 2 Q a {uj), 



N 



1 



y _ 

m\(s + m)! 

m=0 y ' 



sinh 5 

1^ \ s+2m 



(I)' 



(20) 
(21) 

(22) 

(23) 
(24) 



The function £^ d (w;[7) describes the in-scattering of 
polaron-like quasiparticles at the dot (24) . It accounts for 
multiple-phonon emission/absorption processes at finite 
temperature and with finite particle densities. T,^ d (oj; U) 
results from interchanging n B f-> (n B + l), n F -s-> (1— n F ) 
and / O (1 — /) in Eq. (fT5T) . Then the spectral function 
follows using Eqs. (fT2|) and (fl7|) . As we see from Eq. (|T8|) . 
for any 7 < 1 the functions A and / have to be deter- 
mined self-consistently. Moreover, the renormalized dot 
level (|4]) depends on the dot occupation rid, which also 
has to fulfill a self-consistency condition: 



n d 



2^ 



/>; U)A(lu; U). 



(25) 



To determine the variational parameter 7, we minimize 
the thermodynamic potential fl. We use a decoupling ap- 
proximation between the electron and oscillator degrees 
of freedom and neglect the influence of the dot states on 
the leads. We consider an ensemble given by ([3]), but with 
the EP and DL interaction coefficients being multiplied 
by A £ [0, 1]. Then the thermodynamic potential follows 
from the well-known general relations in Rcfs. I22I I25I : 



n 



ln(l + c-^)+e p (l-7)^ 



2^ 



(u-rj)A x {u;U)f x (u;U). (26) 



To make the integration in Eq. (|26j) feasible, we deter- 
mine A\ from Eq. (fTTf with the self-energy in the first 

iteration step, i.e. = A 2 (E^ > + E& }< ). Corre- 

spondingly, f9~' follows from Eq. (p~6|) using S^, and 
r^ 1 ). However, rj will be determined from the dot occu- 
pation rid resulting from the complete self-energy. The 
parameter 7 that minimizes the thermodynamic poten- 
tial determines E^w; U) and, consequently, the com- 
plete functions /(w; U) and A(u; U). 



B. Electron current and linear response 
thermopower 

The operator of the electron current from lead a to the 
dot reads 



3a = 7^ [^^ka 



tka C k a d 



(27) 



with the negative elementary charge e. We determine 
the mean value J a = (J a ) using the connection of the re- 
quired expectation values to the real-time "mixed" Green 
functions g c d{k, a^ti, ti\ U), which are defined similar to 
Eqs. (|6]) and (J7J [21]. In the following we assume identical 
leads and work in the wide band approximation, i.e. we 
set Ta \u>) = Tq = const. Then the steady-state charge 
current through the dot, J = (Jl — Jr)/2, reads 

3 = V / 7T U ) + U L ) - n F (oj + U R )} , 

2 J-00 27T 

(28) 

with the electronic spectral function A(u; U). The latter 
is obtained in terms of the polaronic spectral function as 
follows [H: 



A(u;U) 



(29) 



-7 g coth t 



^I q (k)A{uj;U) + ]T/ S (K)2sinh(s0) 

S>1 

x ( [n B (su ) + /(w + su ; U)] A(u + su ; U) 

+ [n B (sui ) + 1 - /(w - sw ; U)] A(u> - sw ; E/)) j. 

Moreover, we define the differential conductance G of the 
quantum dot system as 



G 



dJ 
d$ 



(30) 



In the linear response regime, we suppose the applica- 
tion of an infinitesimal voltage bias $ = Sfi/e and tem- 
perature difference ST between the leads. Then we can 
expand the current to first order in Sji and ST as [26] 



e 1 



(31) 
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where L is the linear response conductance and X is 
the thermoelectric coefficient. Both quantities follow 
from the linearization of the Fermi functions in Eq. (|28|) 
around the equilibrium chemical potential fj, and the equi- 
librium temperature T: 



L = lim {eJ/Sfi} 

Spi-tO st=o 



X = lim {TJ/5T} 

ST->0 5^=0 



(32) 



eTo 

2 ' 



dw 
2^ 



ui A(u))riF{w)(l — np(ui)). (33) 



In (f3"2")l and (f3"3")l , the electronic spectral function is cal- 
culated in equilibrium. With the help of these transport 
coefficients we define the linear response thermopower 



eX_ 

TL ' 



(34) 



which is a measure of the thermoelectric efficiency of the 
quantum dot system. 



C. Weak EP coupling limit 



The current formula (|28| and the expressions for the 
linear response coefficients in Eqs. ([3"2l and (|33|) have a 
simple structure, because all effects of the EP interac- 
tion are contained in the electronic spectral function A. 
However, our approximation to the spectral function in- 
cludes terms of arbitrarily high order in the EP coupling 
strength g: For 7 > 0, this can be seen explicitly in the 
summations over s in Eqs. (fT9|) and (|29|) . which describe 
inelastic (quasiclastic) processes involving the emission 
and absorption of an unequal (equal) number of phonons. 
As long as 7 < 1, high order terms will also result from 
the iterative calculation of the self-consistent equation 
(|T51) . Via the denominator of the polaronic spectral func- 
tion in Eq. (fT7|) the transport channels will be affected by 
a voltage dependent renormalization of the effective dot 
level and the real part of the self-energy. Lastly, all of 
these contributions are functions of the optimal parame- 
ter 7 m in, which itself will be voltage dependent. This will 
lead to complicated current- voltage characteristics in the 
numerical evaluation of Eq. (|28p , which are presented in 
the next section. 

For a better understanding of the numerical results, we 
want to gain more insight on the different EP coupling 
effects and their dependence on the parameter 7. To this 
end, we consider the limit of small EP coupling strengths 
g and low voltages $ < 2wo- Then we can expand the 
self-energy and the spectral function to second order in 
g around the noninteracting (i.e. zeroth-ordcr) results. 
In doing so, we work in the wide-band approximation 



ri°^(oj) — r = const and consider low temperatures 
T <C u>o, so that «b(wo) ~ 0. First, we set g = in 
Eqs. (|I8p and (|TT)f and obtain the zeroth-order functions 



r(°)( w ) = 2r 0) 



2r 



(35) 
(36) 



(w — A + fi) 2 + Tq 
f \u; U) = \ (n F {u: + U L ) + n F (tu + ET*)) , (37) 



,(0) 



^f(o) {uj] u)A^(u). (38) 



Equations (j 3 5 1) - ([38]) are the exact solution for g = and 
describe a rigid quantum dot acting as a tunneling barrier 
between the leads. Next, we insert A^ and /^°^ for A 
and / in Eq. (fT%|). which corresponds to the first step in 
the self-consistent calculation. Moreover, for T « wo, 
we expand the r.h.s. of Eq. (fl"9|) to second order in g, 
whereby only the terms with s = 0, 1 contribute. The 
resulting approximation of the function T can be written 
as 

r( w; [/)«r(°)H + r( 2 V;io, (39) 

with the second order correction 

r( 2 >( w; u) = -2 7 Vr + 2 7 Vr f/ (0) (w + w o; u) 

+ [(! - 7 )^o] 2 U< >( W + ^o)/ (0) (w + w ; U) 



+A( \u-LUo)(l-f {0) (u-uo;U)) 



(40) 



The second order renormalization of the dot level results 
from substituting n ? for rid in Eq. (|25[) . Then rj is ap- 
proximated as rj rj A — [i + rj^ , with 



n 



(2) =^ p7 (2-7)-2e p (l-7) 2 4 0) 



(41) 



Consequently, we expand the polaronic spectral function 
in Eq. (|17p with respect to the second order corrections 
r^ 2 ' and rj^ and obtain 



with 



A^\u-U) 



A(w) ?s A^(d) + A^(uj; U) , 
' A^\uj) 

2r 



(42) 



x { AT (lu — A + fx) (v {2) + ReE^w; U)) 
+ ((u J -A + rf-T 2 Q y 2 \u J - : U)} (43) 



and 



2tt 



u — w 



(44) 
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Now we replace the polaronic spectral functions on the 
r.h.s. of Eq. with the approximation in Eq. (|4"2"j). 
and keep only terms up to second order in g. Then the 
small coupling approximation to the electronic spectral 
function follows as 

A(u; U) w A<°> (w) - 7 2 .g 2 (w) + A^ 2) (w; 17) 

+ A(°)( w - W o)(l-/ (0) (^-^o;C/))] • (45) 

If we insert T^ 2 ) from Eq. 0DJ into Eq. 03} and sub- 
stitute the resulting expression for A( 2 ) in Eq. (gSJ), the 
electronic spectral function can be written as the sum of 
five terms, 

A(w; U) « A<°>(w) + J4^(w) + 4 2) (w; 17) 

+ 2( 2 )(a;;[/) + 4^(a;;[/), (46) 

whereby A^ ^ is given in Eq. (|36[) and we have defined 



4 2 i» 



7V 



M) 2 



4 2 )(c;[/) = i-(A(°)H) 2 ( W -A + M ) ^ 



(47) 
(48) 



4 2) (o;;C7) 



i-(A(°)H) 2 ( W -A + / i)Re4 2) (c;C/) 



(49) 



+ A(°\u-Lo a )(l-fW(Lo-Lo a ;U)) 



2r n 



(w - A + M ) 2 - r 2 



x {2 7 Vr (/"(°)( W + Wo ;[/) 

+ l-^°)(w-w ;CO) 

+ [(1 - 7 ).9^o] 2 [a<°> ( w + w )/ (0) (w + C/) 

+ A( )( W - Wo )(l-/(°)( W - Wo ; [/))]} . 

(50) 

The function Ap^w) results from the second term on 
the r.h.s. of Eq. (j^J and the first term in Eq. (|3D]|. It 
accounts (to second order) for the polaronic renormal- 
ization of the DL coupling, which gives an overall re- 
duction of the electronic density of states, apart from 

— (21 

the resonance at w = A — fj,. The terms Af,'(u>;U) 

— (2) 

and A^ (w; U) represent the voltage dependent renor- 
malization of the energy levels and contain 7 implicitly. 
Finally, A^ c1 (uj; U) denotes the inelastic contribution to 
the spectral function, which results from tunneling pro- 
cesses that involve the emission of a single phonon at 
the quantum dot. It includes all the terms in the elec- 
tronic spectral function (1451) that contain the functions 



/(°)(w + cj ; U) and 1 - /(°)(w - wo; 17) explicitly. As a 
consequence, it is finite only for |w| > wo and produces 

phononic sidebands in the dot spectrum. However, via 

(2) 

KeT, d J the inelastic channels also contribute to the renor- 

malization of the spectrum at |w| < ujq. Most notably, 
(2) 

for uj —> ±Lu + U a , R- C ^ d causes logarithmic divergences 
in the spectral function. If we evaluate the function /' ' 
for T — s- in Eq. (J40j), then Re£^ 2) follows from Eq. (gl|) 
and contains the logarithmic divergent term 



(1-7) 



2 5 2 _^g£o 

47T 



E 



ln((u-uj + U a ) 2 ) 
(cj - u Q - A + n) 2 + r 2 

ln({uj + uj + U a ) 2 ) \ 
(u + cjo - A + /x) 2 + r 2 J 



(51) 



This term corresponds to the result of Entin-Wohlman et 
al [I?]], but is modified by the prefactor (1 — 7 ) 2 . More- 



over, there is a new contribution to ReS 
term 



(2) 



2t^ 

2 9 1 V"^ 

1 ^T-E 



in 



(OJ - UJQ + Ugf 

(u + Lj + u a y 



namely the 



(52) 



For $ = the logarithmic divergence appearing in 

(21 

ReE^ d (oj; U) for ui — > ojq has the overall prefactor 



4tt 



■(7 2 + 



(l-7) 2 ^ 2 

(A-M) 2 + r 2 



(53) 



so that in the adiabatic (antiadiabatic) limit ujq <C To 
(luq 3> To), an increase in 7 raises (lowers) the overall 
weight of the divergences in the spectral function. 



If we insert Eqs. (|47]) - (f50[) into the current formula 
(|28|) . we get the respective second order corrections to 
the noninteracting current Jw and the differential con- 
ductance, i.e. 



J 
G 



j(o) . 
G (o) 



7(2) 

DL 



G 



(2) 
DL 



4^ 



7(2) 



J; 



(2) 



G 



(2) 



G 



el ' 
(2) 



incl 



(54) 
(55) 



For example, for T — > the second order inelastic tun- 
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ncling current reads 




(56) 



(w + ujo - A + n) 2 + r 2 



+ (l-7)^o 



M) 2 



r 2 
1 



[(^ + W o-A + / i)2 + r2] [( w -A + ^)2 + r2] 



■Ur+uiq 



+ 7^ 



(^-A + M ) 2 -r§ 

[(c-A + M ) 2 +r 2 ] 2 

+ (l-7) 2 ^o 2 



(w - w - a + ii) 2 + r 2 

(w-A + ^-r 2 



[(U-UQ-A + m) 2 + r 2 ] [( w 



A + M ) 2 + r 2 ] 2 



It is finite only for $ > wo , so that the onset of the inelas- 
tic tunneling processes will cause a jump in the differen- 
tial conductance. In general, explicit analytical expres- 
sions for the second order contributions to the differential 
conductance can not be derived, since the optimal param- 
eter 7 m i n is an unknown function of the voltage. However, 
if we suppose that the derivative of 7min ("JO is continuous, 
then for a symmetrical voltage drop Ur — —Ul = e$/2, 
the jump in the differential conductance follows from 
Eq. ([56]) as 



G 



(2) 



incl 



( Up 

\ 1 



a + m) 5 



2tt 



[( 



^A + M ) 2 + r 2 ]- 

2- + (l-7)^ 2 



(57) 



[(-<f -A + / ,) 2 + r 2 ] 

A + m) 2 (-^-A + m ) 2 



\ 2 



r 4 

1 



[(if -A 



M) 2 + r 2 ] 2 [(-f -A + / i) 2 + r 2 ]' 



Again, for 7 — > only the last term on the r.h.s. of 
Eq. (|57|) remains and coincides with the result of Entin- 
Wohlman et al 271. As has been discussed in Ref. ^3, 



this term is negative if the following condition is fulfilled: 



r 2 > 



(A -a*) 5 



(58) 



Then, at $ = ljq, it may cause a downward step in the 
differential conductance. However, the new terms oc 7 2 
in Eq. (|57p are always positive. For large enough 7, they 
outweigh the negative contribution to ([57} . so that the 
overall conductance jumps upwards, even if the condition 
in Eq. §8$ is fulfilled. 



IV. RESULTS AND DISCUSSION 

In the following numerical calculations ujq = 1 is fixed 
as the unit of energy and we set /i = and T = 0.01. 



We work in the wide band approximation, with the 
large bandwidth of the leads W — 60 and T^^w) = 
r 9(u; 2 — (W72) 2 ). The phononic time scale is fixed by 
l/wo, while the electronic time scale is given by 1/ro and 
is used to determine which subsystem is the faster one. 
We will analyze the adiabatic and antiadiabatic limiting 
cases before considering comparable phononic and elec- 
tronic time scales. In doing so, we use the ratio s p /Tq as 
a measure of the EP interaction strength. 

For small to large DL coupling we calculate the po- 
laronic spectral function A and the dot occupation n c i 
self-consistently and determine the variational parame- 
ter 7 ni in by numerically minimizing the thermodynami- 
cal potential 51 as a function of 7. From A, the electronic 
spectral function A as well as the linear response coeffi- 
cients L, X and the particle current J follow. For finite 
voltages, the differential conductance G is calculated nu- 
merically. 

Depending on the bare dot level A, we distinguish be- 
tween the off-resonant (A ^ e p ) and resonant (A = e p ) 
configuration. In the latter case we find that = 0.5 
is a root of (|25|) and we see from Eq. ((4]) that the renor- 
malized dot level resonates with the equilibrium chemical 
potential, i.e. rj = 0, for all 7 m ; n . 



A. Polaron induced NDC 

In their work, La Magna and Deretzis [l5[ suggested 
the variationally determined renormalization of the dot- 
lead coupling as a possible mechanism for the observed 
nonlinear behavior of the differential conductance. We 
investigate whether this remains true within our approx- 
imation, which, in contrast to the effective electron model 
in Ref. Il5l . accounts vibrational features in the electronic 
spectral function to all orders in the EP coupling. 

First we consider the adiabatic regime for weak EP 
coupling by setting To = 10 and e p = 2. We vary the 
voltage bias < $ < 4 and determine the differential 
conductance G. In doing so, we choose an off-resonant 
configuration with A = 8 fixed, so that the dot occu- 
pation is small and remains nearly constant during our 
calculations: m 0.3. 

As a starting point, Fig. [IJa) displays the electronic 
spectral function at <I> = for the variationally deter- 
mined parameter 7 m j n (black line) and compares it to the 
result of a calculation where we kept 7 = fixed instead 
of determining 7 m j n variationally. In general, due to the 
large DL coupling parameter To, the electronic spectral 
function consists of a single wide band. For finite EP 
coupling, vibrational features arise at lo — ±ujq. These 
features can be attributed to logarithmic divergences in 
ReSrf^j as the second order approximation in Eq. (|5ip 
suggests. While they are hardly noticeable for 7 = 0, the 
weight of the logarithmic divergences increases strongly 
in the variational calculation, which yields the optimal 
parameter 7 m in = 0.29. This observation agrees with 
our discussion in the previous section: For the param- 



7 



0.2 



0.15 



~i 1 1 1 1 1 r 



(a) <£=0 



-- y=0 
Y =0.29 

'mm 




FIG. 1: For model parameters T = 0.01, r = 10, e v = 2 and 
A = 8. Panel (a): Electronic spectral functions at $ = for 
fixed 7 = and variationally determined parameter 7 m i n = 
0.29. Panel (b): Differential conductance as a function of 
the voltage bias for 7 = and 7 m m in comparison to the 
noninteracting case e p = 0. Inset: 7 m i n as a function of the 
voltage bias. 
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FIG. 2: For the same parameters as in Fig. [T] The various 
second order contributions to the total differential conduc- 
tance. 



eters used, Eq. (|53|) predicts an increase in the weight 
of the logarithmic contributions by a factor of about 15 
with respect to the 7 = case. Note however, that any 
divergences in the spectrum will be smeared out in our 
results due to the low but finite temperature and a nu- 
merical constraint: We evaluate the self-energy slightly 
above the real uj axis to prevent the unphysical loss of 



spectral weight. 

In Fig. [IJb), the black line presents our result for the 
total differential conductance G as a function of the volt- 
age, with the inset showing the optimal parameter 7 m i n . 
We compare the variational calculation to the cases 7 = 
and e p = 0. For a better understanding of the results 
in Fig. [Ub), the four panels in Fig. [2] show the various 
second order contributions to the total differential con- 
ductance. From Fig. [TJb) it follows that for finite EP 
coupling the overall conductance grows with respect to 
the noninteracting case. Since we are considering the 
off-resonant regime, this can mainly be attributed to the 
lowering of the effective dot level. Accordingly, for 7 = 0, 
we see in Fig. [2] that the function G v accounts for al- 
most all the increase in the conductance. For finite 7 m ; n 

the effective dot level is lowered even further, but the 

(2) 

positive contribution G v is nearly compensated by the 
polaronic renormalization of the DL coupling, which is 
shown in the upper right panel of Fig. [2J With grow- 
ing voltage, the optimal parameter 7 m i n increases. As 
the renormalization of the DL coupling grows stronger, 
a pronounced dip forms in the differential conductance. 
This mechanism is crucial for the interpretation of our 
calculations, as we will see below. 

At <f> = 1, phonon emission by incident electrons be- 
comes possible and opens up an inelastic tunneling chan- 
nel. In the case 7 = 0, we find a small downward step in 
the conductance signal, since with T$ = 10 and A = 8, 
the condition in Eq. ([55)1 is fulfilled. As we discussed in 
the previous section, for finite 7 m i n the first two terms 
on the r.h.s. of Eq. (f57|) can outweigh the third, nega- 
tive term. Accordingly, our numerics show a relatively 
large upward step in the differential conductance (note 
the different scaling factors in the lower right panel of 
Fig. H). 

Next we investigate the polaronic renormalization in 
the antiadiabatic limit (Fq = 0.1) with strong EP cou- 
pling [e p = 2). We choose the resonant configuration 
A = e p . For these parameters, we expect the forma- 
tion of a polaron-like transient state at the quantum dot. 
This is confirmed by the electronic spectral function in 
Fig. Eta), which features several narrow phononic bands. 
In the low-voltage region we find 7, n i n w 0.9, i.e. the 
weight of the variational polaron state is smaller than 
predicted by the complete Lang-Firsov transformation. 

Figure [3]Jb) compares the differential conductance as 
a function of the voltage bias for fixed 7=1 and the 
optimal 7 m j n . Just as in the adiabatic regime consid- 
ered above, we notice small steps in the conductance at 
<!> = 1,3,5 that signal the onset of inelastic transport. 
In addition, a second kind of vibrational feature can be 
found: pronounced conductance peaks arise whenever 
the voltage equals multiple integers of 2wo- Here, res- 
onant transport takes place through the polaronic side 
bands in A. For 7 = 1 the differential conductance 
stays strictly positive, but approaches zero between these 
well separated peaks. As seen for the adiabatic case, in 
the full calculation the polaronic renormalization grows 
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FIG. 3: For model parameters T = 0.01, To = 0.1, e v = 2 and 
A = 2. Panel (a): Electronic spectral function of the varia- 
tional calculation for $ = and 7 m m = 0.9. Panel (b): Differ- 
ential conductance as a function of the voltage bias, compared 
to the result with fixed 7 = 1. 
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FIG. 4: For model parameters T = 0.01, r = 1, A = 2 
and e p = 2. Panel (a): Electronic spectral function in the 
variational calculation for $ = and 7 m in = 0.5. Panel (b): 
Differential conductance as a function of the voltage for the 
variational calculation (7 m i n ), compared to the results with 
fixed 7 = and 7 = 1. Inset: Optimized variational parame- 
ter as a function of the voltage. 



stronger with increasing voltage bias. As a consequence, 
in the low voltage region the differential conductance be- 
comes negative between the resonance peaks. Note how- 
ever, that at $ = 1 and $ = 3 the positive nonresonant 
conductance steps, although carrying little weight, ren- 
der the differential conductance positive again. 

Thanks to our variational approach, we are able to in- 
vestigate the interesting regime of comparable electronic 
and phonemic energies. To this end, we set To = 1 and 
consider intermediate EP coupling e v = 2. As before, 
we examine the resonant, electron-holc-symmetric situa- 
tion with A = 2. Fig@Ja) shows the electronic spectral 
function at zero voltage, where the variational calcula- 
tion yields 7 m i n ~ 0.5. Due to comparable electronic and 
phononic time scales, the width of the few phononic side 
bands is of the order of their spacing. 

In Fig. |U[b) we compare the conductance signal of the 
variational calculation to both, the 7 = and 7 = 1 
cases. In the low voltage regime, we have 7 m i n > 0.5 
and the DL coupling is moderately rcnormalizcd. As 
the voltage grows, the variational parameter steadily in- 
creases and, as can be seen in the inset of Fig@Jb), the 
polaron effect strengthens whenever a new resonant in- 
elastic channel is accessible. The vibrational features in 
the conductance signal are heavily modulated by the volt- 
age dependent polaronic renormalization: In contrast to 
the cases with fixed 7, there is no clear distinction be- 
tween resonant peaks at $ = 2,4 and off- resonant steps 
at $ = 1,3,5, since the latter become peaks, too. Due 
to the comparable phononic and electronic time scales, 



both kinds of vibrational features have nearly the same 
spectral weight. Moreover, the differential conductance 
approaches zero between the broad conductance peaks, 
but no NDC is observed. 

To sum up, the polaron formation involves the redis- 
tribution of spectral weight in the local density of states 
and, most importantly, the renormalization of the effec- 
tive DL coupling. For strong EP interaction, it is indeed a 
possible mechanism for NDC. Yet, for small to intermedi- 
ate coupling, the NDC is suppressed when multi-phonon 
transport processes are taken into account. 

B. Effective dot level 

In the following we present another interesting con- 
sequence of the variational polaronic renormalization, 
which concerns the effective dot level. 

We choose a slightly off-resonant configuration with 
A = and e p > 0, so that in contrast to the above 
calculations, the effective dot level is not pinned to the 
equilibrium chemical potential. We decrease the bare DL 
coupling slightly (Fo = 0.33) and consider weak to inter- 
mediate EP coupling strengths. Figure [SJa) compares 
the differential conductance for 7 = to the variational 
calculation. For weak EP coupling, Fig. E{b) shows the 

second order approximations and G^L. Fig. [3Jc) 
finally presents the optimal parameter 7 m j n and the ef- 
fective dot level rj as functions of the voltage. 
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FIG. 5: For T = 0.01, r = 0.33 and A = 0. Panel (a): 
Differential conductance as a function of the voltage for the 
variational calculation (7 m i n ), compared to the results with 
fixed 7 = 0. Panel (b): Second order contributions to the dif- 
ferential conductance for e p — 0.11. Panel (c): Variationally 
determined parameter 7 m i n and renormalized dot-level rj as 
functions of the voltage for e p = 0.33. 



With e p = 0.11, the system parameters correspond to 
the case of high zeroth order transmission presented in 
Fig. 5(a) in the work of Entin-Wohlman et al 12711. Our 
result for 7 = is in good agreement with Ref. [27J. The 
conductance maximum lies near $ = and we hnd a 
small conductance dip at $ = 1, which is caused by the 
logarithmic divergence in RcS^- In the full calculation 
for e p = 0.11, we find 7 m j n = 0.75 at $ = 1. Here 
the dip in the total conductance vanishes. The second 
order approximation in the left panel of Fig.[SJb) suggests 
that this is mainly due to a reduction of the weight of 
the logarithmic divergence in ReS^d- From Fig. [5jb) we 

(2) 

can also see that the jump in G incl is positive for both 
7 = and 7 m i n , since the condition (|58p is not fulfilled 
for the given parameters. Moreover, in the variational 
calculation the height of the conductance jump is reduced 



with respect to the 7 = result, which can be confirmed 
using Eq. ([571). 

If we increase the EP coupling to e v = 0.33, the dip 
in the total conductance reappears. But most impor- 
tantly, instead of a broad conductance resonance, we 
find a peak-like feature at $ = 0.7. As we see from 
Fig. (njc), with increasing voltage rj shifts upwards until 
at $ = 0.62 it approaches the chemical potential of a 
lead. For 0.65 < $ < 0.8, the variational parameter de- 
creases in such a way that rj stays in resonance with the 
lead chemical potential. The decrease in 7 m i n reduces the 
renormalization of the DL coupling. Thereby, the system 
maximizes the resonant tunneling current with respect to 
the 7 = case and a new peak-like conductance feature 
is observed in Fig. [5{a). This "sticking" of the effective 
dot level to the lead chemical potentials is the second 
main result of our variational calculations. 

Now we consider the off-resonant scenario A = 10 for 
strong EP coupling e p = 8. The results are presented 
in Fig. |5] (note that Fig. [SJa) shows the total current). 
As expected, the effective dot level rj sinks notably with 
growing voltage, until at $ = 6.2 it begins to grow lin- 
early, following the upper lead chemical potential. Again, 
the differential conductance grows considerably. In con- 
trast to the intermediate EP coupling case, 7 m i„ jumps 
from 0.4 to 0.6 when the system switches between two 
local minima in the thermodynamic potential. The re- 
sulting discontinuities in rj and Tq cause a noticeable drop 
in the total current. As the voltage grows further, 7 m i n 
decreases again. Now the first phonon side band at rj+ujQ 
sticks to the lead chemical potential and the conductance 
grows once more. Similar behavior, involving the second 
and third side bands, is found at 4> ~ 9 and <f> ps 10.5, re- 
spectively, until 7 m in = 1 in the high- voltage limit. More- 
over, due to the strong EP coupling, the upward steps in 
the current are followed by regions with NDC. 



C. Thermopower 

Finally, we investigate the thermoelectric response of 
the molecular junction in the physically most interesting 
regime of intermediate DL coupling. Setting To = 1, we 
consider the equilibrium situation $ = 0. For e p = 2, 
we compare the variational calculation to the cases with 
fixed 7 = 0,1 and to the noninteracting system e p = 0. 
Fig. [TJa) shows the linear response thermopower S as a 
function of the bare dot level, while Fig.[7][b) presents the 
thermoelectric coefficient and the linear response conduc- 
tance. 

In general, S features two resonances of opposite sign. 
For £ p = 2 they are located at A ss e v ± 0.2. In the small 
polaron limit 7=1 our calculation predicts a substantial 
increase in the maximum thermopower with respect to 
all the other cases. This can be explained with the help 
of the respective electronic spectral functions plotted in 
Fig. 0(c) for A = 2.2. In the case 7 = 0, the spectral 
function features a broad band around the Fermi edge at 
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FIG. 6: For model parameters T = 0.01, r = 1, A = 10 
and e p = 8. Panel (a): Electron current as a function of the 
voltage. Panel (b): Variationally determined parameter 7 m i n 
and renormalized dot-level rj as functions of the voltage. 



io = 0. The states with high energies uj > have only 
slightly more spectral weight than the states with low 
energies ui < 0. Because the integrand on the r.h.s. of 
Eq. ([55)1 is weighted by u>, the resulting thermoelectric 
response coefficient X is small. Physically, this means 
that a small temperature difference between the leads in- 
duces the flow of high energy particles through the quan- 
tum dot, which, in principle, can result in a voltage drop 
across the junction. In the case 7 = however, the cur- 
rent is compensated by a nearly equal counterflow of low 
energy carriers, so that the overall thermoelectric effect 
is small. If A is lowered to 1.8, the low-energy states 
have the larger spectral weight and the thermoelectric 
response coefficient changes sign. For A = e p = 2, the 
spectral function is symmetric around ui = so that the 
net charge current induced by the temperature difference 
vanishes and we have X = 0. 

For 7 = 1, the strong renormalization of the DL cou- 
pling reduces the width of the bands in the spectral func- 
tion in Fig. [TJc). As a result, near the Fermi edge the 
relative weight of the high-energy states increases, so that 
the dot acts as a more effective energy filter. The unfa- 
vorable counterflow of low-energy charge carriers is sup- 
pressed and the thermoelectric response X grows consid- 
erably (see Fig. J7](b)). As can also be seen in Fig. [TJb), 
the linear response conductance L in Eq. (|52"|) decreases 
when 7 is set from zero to one, since it depends only on 




FIG. 7: For model parameters T = 0.01, r = 1, e p = 2 
and $ = 0. Panel (a): Thermopower as a function of the 
bare dot level in the noninteracting (e p = 0) and interacting 
system. Panel (b): Thermoelectric response X and linear 
conductance L as functions of the bare dot level. Panel (c): 
Electronic spectral functions at A = 2.2 for 7 = 0, 1 and 

7min = 0.5. 



the (shrinking) spectral weight around the Fermi edge. 
This, too, boosts the thermopower S. 

At A = 2.2 the variational calculation yields 7 m i n = 
0.5, so that the width of the zcro-phonon band lies be- 
tween the other results. Consequently, this is also true for 
the maximum value of X. Note however, that our varia- 
tional calculation maximizes the linear response conduc- 
tance L with respect to both limiting cases, so that the 
maximum thermopower is only slightly larger than for 
7 = 0. We conclude, that the local EP interaction can, 
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in principle, enhance the maximum thcrmopower of the 
quantum dot device. Yet, for intermediate DL coupling 
strengths the small polaron picture with 7 = 1 greatly 
overestimates the effect. 



V. CONCLUDING REMARKS 

To summarize, adopting a generalized variational 
Lang-Firsov transformation we calculate the interacting 
spectral function of a molecular quantum dot for small- 
to-largc DL coupling and weak-to-strong EP interaction. 
We investigate the impact of the formation of a polaronic 
dot state on the steady-state current-voltage character- 
istics, as well as on the linear response thcrmopower of 
the system. 

In the case of strong EP interaction, the voltage- 
dependent polaronic renormalization of the DL coupling 
causes negative differential conductance. For compara- 
ble electronic and phononic time scales, this effect is di- 
minished by transport through overlapping phonon side 
bands. 

We find that in the off-resonant or ungated configura- 



tion, the renormalized dot level follows the lead chemical 
potentials. This process generates new peaks in the dif- 
ferential conductance signal. 

In the equilibrium situation, the EP coupling enhances 
the thermopower of the quantum dot device, albeit by a 
smaller factor than predicted in the small polaron limit. 

The present work may be extended in several direc- 
tions. Most notably, in the nonequilibrium regime, one 
should investigate the impact of the observed NDC on the 
thermoelectric properties of the molecular junction. The 
dynamics and heating of the vibrational subsystem could 
be included by means of nonequilibrium phonon Green 
functions (28|. Moreover, in the light of recent advances 
in nanotechnology and experimental studies, new geome- 
tries have come into focus, like multi-terminal junctions 
or a molecule placed on an Aharonov-Bohm ring [29| . 
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